# Perceived Discrimination and Political Behavior Replication Code
# British Journal of Political Science
# Date: Feb. 6, 2018
# Author: Kassra A.R. Oskooii (oskooiik@udel.edu)


rm(list=ls())

#Load Packages (Install Needed Packages)
library(MASS)
library(simcf) #Download Package: https://faculty.washington.edu/cadolph/index.php?page=60
library(tile) #Download Package: https://faculty.washington.edu/cadolph/index.php?page=60
library(nlme)
library(mgcv)
library(RColorBrewer)
library(stargazer)
library(psych)
library(dplyr)


#Load Dataset
data<-haven::read_dta("EMBES_BJPS_DATA.dta")


######################
##Tables 1 & 2 Set-up
######################
###################################
##GENERAL ELECTION (Broad Disc I)
###################################

formula<-(voted2010 ~ 
            disc_race_binary + 
          relatt_oth_r + pol_interest + polknowledge + partyid +
          identity+ english+ native_born+ female+ age+ education2+ highinc+ 
          medinc+ misinc+ black_caribbean+ indian+ pakistani+ bangladeshi+
          vote_duty+ efficacy+ democ_satis+ trust_parliament)

mdata<-extractdata(formula, data, na.rm=TRUE)

k <- length(names(mdata))
for (k in 1:k) {
  print(class(mdata[,k]))
}
for (k in 1:k) {
  mdata[,k] <- as.numeric(mdata[,k])
}

logit.result1<- glm(formula, family=binomial(link=logit), data=mdata)
summary(logit.result1)

pe.1 <- logit.result1$coefficients  # point estimates
vc.1 <- vcov(logit.result1)

### simulation
sims <- 10000
simbetas <- mvrnorm(sims, pe.1, vc.1)

xhyp1 <- cfMake(formula, mdata, nscen=1, f=mean)

xhyp1 <- cfName(xhyp1, "disc_race_binary", scen=1) 
xhyp1 <-cfChange(xhyp1, "disc_race_binary", x=1, xpre=0, scen=1)

xhyp1
yhyp1<-logitsimfd(xhyp1, simbetas, ci=0.90)
yhyp1


###################################
###LOCAL Election (Broad Disc I)
###################################

formula<-(voted2010_local ~ 
            disc_race_binary + 
            relatt_oth_r + pol_interest + polknowledge + partyid +
            identity+ english+ native_born+ female+ age+ education2+ highinc+ 
            medinc+ misinc+ black_caribbean+ indian+ pakistani+ bangladeshi+
            vote_duty+ efficacy+ democ_satis+ trust_parliament)


mdata<-extractdata(formula, data, na.rm=TRUE)

k <- length(names(mdata))
for (k in 1:k) {
  print(class(mdata[,k]))
}
for (k in 1:k) {
  mdata[,k] <- as.numeric(mdata[,k])
}

logit.result2<- glm(formula, family=binomial(link=logit), data=mdata)
summary(logit.result2)

pe.2 <- logit.result2$coefficients  # point estimates
vc.2 <- vcov(logit.result2)

### simulation
sims <- 10000
simbetas <- mvrnorm(sims, pe.2, vc.2)

xhyp2 <- cfMake(formula, mdata, nscen=1, f=mean)

xhyp2 <- cfName(xhyp2, "disc_race_binary", scen=1) 
xhyp2 <-cfChange(xhyp2, "disc_race_binary", x=1, xpre=0, scen=1)

xhyp2
yhyp2<-logitsimfd(xhyp2, simbetas, ci=0.90)
yhyp2



###################################
###GENERAL Election (Broad Disc II)
###################################

formula<-(voted2010 ~ 
            disc_race +
            relatt_oth_r + pol_interest + polknowledge + partyid +
            identity+ english+ native_born+ female+ age+ education2+ highinc+ 
            medinc+ misinc+ black_caribbean+ indian+ pakistani+ bangladeshi+
            vote_duty+ efficacy+ democ_satis+ trust_parliament)


mdata<-extractdata(formula, data, na.rm=TRUE)

k <- length(names(mdata))
for (k in 1:k) {
  print(class(mdata[,k]))
}
for (k in 1:k) {
  mdata[,k] <- as.numeric(mdata[,k])
}

logit.result1.1<- glm(formula, family=binomial(link=logit), data=mdata)
summary(logit.result1.1)

pe.1.1 <- logit.result1.1$coefficients  # point estimates
vc.1.1 <- vcov(logit.result1.1)

### simulation
sims <- 10000
simbetas <- mvrnorm(sims, pe.1.1, vc.1.1)

xhyp1.1 <- cfMake(formula, mdata, nscen=1, f=mean)

xhyp1.1 <- cfName(xhyp1.1, "disc_race", scen=1) 
xhyp1.1 <-cfChange(xhyp1.1, "disc_race", x=3, xpre=0, scen=1)

xhyp1.1
yhyp1.1<-logitsimfd(xhyp1.1, simbetas, ci=0.90)
yhyp1.1



###################################
###LOCAL Election (Broad Disc II)
###################################

formula<-(voted2010_local ~ 
            disc_race + 
            relatt_oth_r + pol_interest + polknowledge + partyid +
            identity+ english+ native_born+ female+ age+ education2+ highinc+ 
            medinc+ misinc+ black_caribbean+ indian+ pakistani+ bangladeshi+
            vote_duty+ efficacy+ democ_satis+ trust_parliament)


mdata<-extractdata(formula, data, na.rm=TRUE)

k <- length(names(mdata))
for (k in 1:k) {
  print(class(mdata[,k]))
}
for (k in 1:k) {
  mdata[,k] <- as.numeric(mdata[,k])
}

logit.result2.1<- glm(formula, family=binomial(link=logit), data=mdata)
summary(logit.result2.1)

pe.2.1 <- logit.result2.1$coefficients  # point estimates
vc.2.1 <- vcov(logit.result2.1)

### simulation
sims <- 10000
simbetas <- mvrnorm(sims, pe.2.1, vc.2.1)

xhyp2.1 <- cfMake(formula, mdata, nscen=1, f=mean)

xhyp2.1 <- cfName(xhyp2.1, "disc_race", scen=1) 
xhyp2.1 <-cfChange(xhyp2.1, "disc_race", x=3, xpre=0, scen=1)

xhyp2.1
yhyp2.1<-logitsimfd(xhyp2.1, simbetas, ci=0.90)
yhyp2.1



######################################
#GENERAL Election (Specific Measures)
######################################

formula<-(voted2010 ~ 
            socdisc + poldisc+
            relatt_oth_r + pol_interest + polknowledge + partyid +
            identity+ english+ native_born+ female+ age+ education2+ highinc+ 
            medinc+ misinc+ black_caribbean+ indian+ pakistani+ bangladeshi+
            vote_duty+ efficacy+ democ_satis+ trust_parliament)

mdata<-extractdata(formula, data, na.rm=TRUE)

k <- length(names(mdata))
for (k in 1:k) {
  print(class(mdata[,k]))
}
for (k in 1:k) {
  mdata[,k] <- as.numeric(mdata[,k])
}

logit.result1.2<- glm(formula, family=binomial(link=logit), data=mdata)
summary(logit.result1.2)

pe.1.2 <- logit.result1.2$coefficients  # point estimates
vc.1.2 <- vcov(logit.result1.2)

### simulation
sims <- 10000
simbetas <- mvrnorm(sims, pe.1.2, vc.1.2)

xhyp1.2 <- cfMake(formula, mdata, nscen=2, f=mean)

xhyp1.2 <- cfName(xhyp1.2, "socdisc", scen=1) 
xhyp1.2 <-cfChange(xhyp1.2, "socdisc", x=9, xpre=0, scen=1)

xhyp1.2 <- cfName(xhyp1.2, "poldisc", scen=2) 
xhyp1.2 <-cfChange(xhyp1.2, "poldisc", x=12, xpre=0, scen=2)

xhyp1.2
yhyp1.2<-logitsimfd(xhyp1.2, simbetas, ci=0.90)
yhyp1.2



####################################
#LOCAL Election (Specific Measures)
####################################

formula<-(voted2010_local ~ 
            socdisc + poldisc+
            relatt_oth_r + pol_interest + polknowledge + partyid +
            identity+ english+ native_born+ female+ age+ education2+ highinc+ 
            medinc+ misinc+ black_caribbean+ indian+ pakistani+ bangladeshi+
            vote_duty+ efficacy+ democ_satis+ trust_parliament)


mdata<-extractdata(formula, data, na.rm=TRUE)

k <- length(names(mdata))
for (k in 1:k) {
  print(class(mdata[,k]))
}
for (k in 1:k) {
  mdata[,k] <- as.numeric(mdata[,k])
}

logit.result2.2<- glm(formula, family=binomial(link=logit), data=mdata)
summary(logit.result2.2)

pe.2.2 <- logit.result2.2$coefficients  # point estimates
vc.2.2 <- vcov(logit.result2.2)

### simulation
sims <- 10000
simbetas <- mvrnorm(sims, pe.2.2, vc.2.2)

xhyp2.2 <- cfMake(formula, mdata, nscen=2, f=mean)

xhyp2.2 <- cfName(xhyp2.2, "socdisc", scen=1) 
xhyp2.2 <-cfChange(xhyp2.2, "socdisc", x=9, xpre=0, scen=1)

xhyp2.2 <- cfName(xhyp2.2, "poldisc", scen=2) 
xhyp2.2 <-cfChange(xhyp2.2, "poldisc", x=12, xpre=0, scen=2)

xhyp2.2
yhyp2.2<-logitsimfd(xhyp2.2, simbetas, ci=0.90)
yhyp2.2



######################################
#GENERAL Election (Simplified Model)
######################################

formula<-(voted2010 ~ 
            socdisc + poldisc+
            relatt_oth_r + pol_interest + polknowledge + partyid +
            identity+ english+ native_born+ female+ age+ education2+ highinc+ 
            medinc+ misinc+ black_caribbean+ indian+ pakistani+ bangladeshi)


mdata<-extractdata(formula, data, na.rm=TRUE)

k <- length(names(mdata))
for (k in 1:k) {
  print(class(mdata[,k]))
}
for (k in 1:k) {
  mdata[,k] <- as.numeric(mdata[,k])
}

result_simple_general<- glm(formula, family=binomial(link=logit), data=mdata)
summary(result_simple_general)



######################################
#LOCAL Election (Simplified Model)
######################################

formula<-(voted2010_local ~ 
            socdisc + poldisc+
            relatt_oth_r + pol_interest + polknowledge + partyid +
            identity+ english+ native_born+ female+ age+ education2+ highinc+ 
            medinc+ misinc+ black_caribbean+ indian+ pakistani+ bangladeshi)


mdata<-extractdata(formula, data, na.rm=TRUE)

k <- length(names(mdata))
for (k in 1:k) {
  print(class(mdata[,k]))
}
for (k in 1:k) {
  mdata[,k] <- as.numeric(mdata[,k])
}

result_simple_local<- glm(formula, family=binomial(link=logit), data=mdata)
summary(result_simple_local)


###################
## Variable Lables 
###################
modellabels<-c("Borad Disc I", "Broad Disc II", "Societal Disc", "Political Disc",
               "Worship Attendance", "Political Interest", "Political Knowledge", "Party ID (Yes=1)", 
               "Identity (Brit=2)", "English (Main Lang)", "Native Born", "Female", 
               "Age", "Education", "High Income", "Med Income", "Missing Income",
               "Black Caribbean", "Indian", "Pakistani", "Bangladeshi",
               "Vote Duty", "Political Efficacy","Democratic Satisfaction","Trust Parliament")


dvlabel1<-c("Vote in General Election")
dvlabel2<-c("Vote in Local Election")


##########
##TABLE 1
##########
stargazer(logit.result2, logit.result2.1, logit.result2.2, result_simple_local, style="APSR", covariate.labels = modellabels, out.header=T,
          model.numbers = TRUE, 
          title = "",
          dep.var.labels = dvlabel2)

##########
##TABLE 2
##########
stargazer(logit.result1, logit.result1.1, logit.result1.2, result_simple_general, style="APSR", covariate.labels = modellabels, out.header=T,
          model.numbers = TRUE, 
          title = "",
          dep.var.labels = dvlabel1)


###########################################
####Figure 3: Local Election (LEFT PANEL)
###########################################

#Local Election 1
trace1 <- ropeladder(
  x=yhyp2$pe,
  lower=yhyp2$lower,
  upper=yhyp2$upper,
  labels=c("Disc. (0-1)"), 
  col= c("grey40"),
  pch=c(16),
  cex=2, 
  plot=1)

#Local Election 2
trace2 <- ropeladder(
  x=yhyp2.1$pe,
  lower=yhyp2.1$lower,
  upper=yhyp2.1$upper,
  labels=c("Disc. (0-3)"), 
  col= c("grey40"),
  pch=c(16),
  cex=2, 
  plot=2)

#Local Election 3
trace3 <- ropeladder(
  x=yhyp2.2$pe,
  lower=yhyp2.2$lower,
  upper=yhyp2.2$upper,
  labels=c("Societal Disc.", "Political Disc."), 
  col= c("Black", "Black"),
  pch=c(16, 16),
  cex=2, 
  plot=3)

vertmark <- linesTile(x = c(0,0),
                      y = c(0,1),
                      lty = "dashed",
                      col= "Black",
                      plot = 1:3
)

# Widen space between titles on left
trace1$entryheight <- .25
trace2$entryheight <- .25
trace3$entryheight <- .15

at.x1 <- seq(-.2, .2,.1)
at.x2 <- seq(-.2, .2,.1)
at.x3 <- seq(-.4, .4,.1)

tc<-tile(trace1,
         trace2,
         trace3,
         vertmark,
         RxC=c(3,1),
         xaxis=list(at1 = at.x1, 
                    at2= at.x2, 
                    at3= at.x3, 
                    cex=.8),
         topaxis=list(at1 = at.x1, 
                      at2= at.x2, 
                      at3= at.x3,  
                      cex=.8,
                      add = rep(TRUE,3)),
         topaxistitle = list(labels=c("Local Election", "", ""), cex=1),
         output = list(file = "EMBES_figure3_local", type="pdf", width=10), 
         frame=FALSE,
         gridlines=list(type="t", col="grey65")
)



##############################################
####Figure 3: General Election (RIGHT PANEL)
##############################################

#General Election 1
trace1 <- ropeladder(
  x=yhyp1$pe,
  lower=yhyp1$lower,
  upper=yhyp1$upper,
  labels=c("Disc. (0-1)"), 
  col= c("grey40"),
  pch=c(16),
  cex=2, 
  plot=1)

#General Election 2
trace2 <- ropeladder(
  x=yhyp1.1$pe,
  lower=yhyp1.1$lower,
  upper=yhyp1.1$upper,
  labels=c("Disc. (0-3)"), 
  col= c("grey40"),
  pch=c(16),
  cex=2, 
  plot=2)

#General Election 3
trace3 <- ropeladder(
  x=yhyp1.2$pe,
  lower=yhyp1.2$lower,
  upper=yhyp1.2$upper,
  labels=c("Societal Disc.", "Political Disc."), 
  col= c("Black", "Black"),
  pch=c(16, 16),
  cex=2, 
  plot=3)

vertmark <- linesTile(x = c(0,0),
                      y = c(0,1),
                      lty = "dashed",
                      col= "Black",
                      plot = 1:3
)

# Widen space between titles on left
trace1$entryheight <- .25
trace2$entryheight <- .25
trace3$entryheight <- .15

at.x1 <- seq(-.2, .2,.1)
at.x2 <- seq(-.2, .2,.1)
at.x3 <- seq(-.4, .4,.1)


tc<-tile(trace1,
         trace2,
         trace3,
         vertmark,
         RxC=c(3,1),
         xaxis=list(at1 = at.x1, 
                    at2= at.x2, 
                    at3= at.x3, 
                    cex=.8),
         topaxis=list(at1 = at.x1, 
                      at2= at.x2, 
                      at3= at.x3,  
                      cex=.8,
                      add = rep(TRUE,3)),
         topaxistitle = list(labels=c("General Election", " ", " "), cex=1),
         output = list(file = "EMBES_figure3_general", type="pdf", width=10), 
         frame=FALSE,
         gridlines=list(type="t", col="grey65")
)




#####################################################################
##Figure 4: Predictors of Voting in the Local and General Elections
#####################################################################

#Local
formula<-(voted2010_local ~ 
            socdisc + poldisc+
            relatt_oth_r + pol_interest + polknowledge + partyid +
            identity+ english+ native_born+ female+ age+ education2+ highinc+ 
            medinc+ misinc+ black_caribbean+ indian+ pakistani+ bangladeshi+
            vote_duty+ efficacy+ democ_satis+ trust_parliament)


mdata<-extractdata(formula, data, na.rm=TRUE)

k <- length(names(mdata))
for (k in 1:k) {
  print(class(mdata[,k]))
}
for (k in 1:k) {
  mdata[,k] <- as.numeric(mdata[,k])
}

logit.result2.2<- glm(formula, family=binomial(link=logit), data=mdata)
summary(logit.result2.2)

pe.2.2 <- logit.result2.2$coefficients  
vc.2.2 <- vcov(logit.result2.2)

# simulation
sims <- 10000
simbetas <- mvrnorm(sims, pe.2.2, vc.2.2)

xhyp2.2 <- cfMake(formula, mdata, nscen=22, f=mean)

xhyp2.2 <- cfName(xhyp2.2, "socdisc", scen=1) 
xhyp2.2 <-cfChange(xhyp2.2, "socdisc", x=9, xpre=0, scen=1)

xhyp2.2 <- cfName(xhyp2.2, "poldisc", scen=2) 
xhyp2.2 <-cfChange(xhyp2.2, "poldisc", x=12, xpre=0, scen=2)

xhyp2.2 <- cfName(xhyp2.2, "relatt_oth_r", scen=3) 
xhyp2.2 <-cfChange(xhyp2.2, "relatt_oth_r", x=5, xpre=0, scen=3)

xhyp2.2 <- cfName(xhyp2.2, "pol_interest", scen=4) 
xhyp2.2 <-cfChange(xhyp2.2, "pol_interest", x=4, xpre=0, scen=4)

xhyp2.2 <- cfName(xhyp2.2, "polknowledge", scen=5) 
xhyp2.2 <-cfChange(xhyp2.2, "polknowledge", x=5, xpre=0, scen=5)

xhyp2.2 <- cfName(xhyp2.2, "partyid", scen=6) 
xhyp2.2 <-cfChange(xhyp2.2, "partyid", x=1, xpre=0, scen=6)

xhyp2.2 <- cfName(xhyp2.2, "identity", scen=7) 
xhyp2.2 <-cfChange(xhyp2.2, "identity", x=2, xpre=0, scen=7)

xhyp2.2 <- cfName(xhyp2.2, "english", scen=8) 
xhyp2.2 <-cfChange(xhyp2.2, "english", x=1, xpre=0, scen=8)

xhyp2.2 <- cfName(xhyp2.2, "native_born", scen=9) 
xhyp2.2 <-cfChange(xhyp2.2, "native_born", x=1, xpre=0, scen=9)

xhyp2.2 <- cfName(xhyp2.2, "female", scen=10) 
xhyp2.2 <-cfChange(xhyp2.2, "female", x=1, xpre=0, scen=10)

xhyp2.2 <- cfName(xhyp2.2, "age", scen=11) 
xhyp2.2 <-cfChange(xhyp2.2, "age", x=40, xpre=0, scen=11)

xhyp2.2 <- cfName(xhyp2.2, "education2", scen=12) 
xhyp2.2 <-cfChange(xhyp2.2, "education2", x=5, xpre=0, scen=12)

xhyp2.2 <- cfName(xhyp2.2, "highinc", scen=13) 
xhyp2.2 <-cfChange(xhyp2.2, "highinc", x=1, xpre=0, scen=13)

xhyp2.2 <- cfName(xhyp2.2, "medinc", scen=14) 
xhyp2.2 <-cfChange(xhyp2.2, "medinc", x=1, xpre=0, scen=14)

xhyp2.2 <- cfName(xhyp2.2, "black_caribbean", scen=15) 
xhyp2.2 <-cfChange(xhyp2.2, "black_caribbean", x=1, xpre=0, scen=15)

xhyp2.2 <- cfName(xhyp2.2, "indian", scen=16) 
xhyp2.2 <-cfChange(xhyp2.2, "indian", x=1, xpre=0, scen=16)

xhyp2.2 <- cfName(xhyp2.2, "pakistani", scen=17) 
xhyp2.2 <-cfChange(xhyp2.2, "pakistani", x=1, xpre=0, scen=17)

xhyp2.2 <- cfName(xhyp2.2, "bangladeshi", scen=18) 
xhyp2.2 <-cfChange(xhyp2.2, "bangladeshi", x=1, xpre=0, scen=18)

xhyp2.2 <- cfName(xhyp2.2, "vote_duty", scen=19) 
xhyp2.2 <-cfChange(xhyp2.2, "vote_duty", x=5, xpre=1, scen=19)

xhyp2.2 <- cfName(xhyp2.2, "efficacy", scen=20) 
xhyp2.2 <-cfChange(xhyp2.2, "efficacy", x=10, xpre=0, scen=20)

xhyp2.2 <- cfName(xhyp2.2, "democ_satis", scen=21) 
xhyp2.2 <-cfChange(xhyp2.2, "democ_satis", x=3, xpre=0, scen=21)

xhyp2.2 <- cfName(xhyp2.2, "trust_parliament", scen=22) 
xhyp2.2 <-cfChange(xhyp2.2, "trust_parliament", x=10, xpre=0, scen=22)

xhyp2.2
yhyp2.2<-logitsimfd(xhyp2.2, simbetas, ci=0.90)
yhyp2.2


#General
formula<-(voted2010 ~ 
            socdisc + poldisc+
            relatt_oth_r + pol_interest + polknowledge + partyid +
            identity+ english+ native_born+ female+ age+ education2+ highinc+ 
            medinc+ misinc+ black_caribbean+ indian+ pakistani+ bangladeshi+
            vote_duty+ efficacy+ democ_satis+ trust_parliament)

mdata<-extractdata(formula, data, na.rm=TRUE)

k <- length(names(mdata))
for (k in 1:k) {
  print(class(mdata[,k]))
}
for (k in 1:k) {
  mdata[,k] <- as.numeric(mdata[,k])
}

logit.result1.2<- glm(formula, family=binomial(link=logit), data=mdata)
summary(logit.result1.2)

pe.1.2 <- logit.result1.2$coefficients 
vc.1.2 <- vcov(logit.result1.2)

# simulation
sims <- 10000
simbetas <- mvrnorm(sims, pe.1.2, vc.1.2)

xhyp1.2 <- cfMake(formula, mdata, nscen=22, f=mean)

xhyp1.2 <- cfName(xhyp1.2, "socdisc", scen=1) 
xhyp1.2 <-cfChange(xhyp1.2, "socdisc", x=9, xpre=0, scen=1)

xhyp1.2 <- cfName(xhyp1.2, "poldisc", scen=2) 
xhyp1.2 <-cfChange(xhyp1.2, "poldisc", x=12, xpre=0, scen=2)

xhyp1.2 <- cfName(xhyp1.2, "relatt_oth_r", scen=3) 
xhyp1.2 <-cfChange(xhyp1.2, "relatt_oth_r", x=5, xpre=0, scen=3)

xhyp1.2 <- cfName(xhyp1.2, "pol_interest", scen=4) 
xhyp1.2 <-cfChange(xhyp1.2, "pol_interest", x=4, xpre=0, scen=4)

xhyp1.2 <- cfName(xhyp1.2, "polknowledge", scen=5) 
xhyp1.2 <-cfChange(xhyp1.2, "polknowledge", x=5, xpre=0, scen=5)

xhyp1.2 <- cfName(xhyp1.2, "partyid", scen=6) 
xhyp1.2 <-cfChange(xhyp1.2, "partyid", x=1, xpre=0, scen=6)

xhyp1.2 <- cfName(xhyp1.2, "identity", scen=7) 
xhyp1.2 <-cfChange(xhyp1.2, "identity", x=2, xpre=0, scen=7)

xhyp1.2 <- cfName(xhyp1.2, "english", scen=8) 
xhyp1.2 <-cfChange(xhyp1.2, "english", x=1, xpre=0, scen=8)

xhyp1.2 <- cfName(xhyp1.2, "native_born", scen=9) 
xhyp1.2 <-cfChange(xhyp1.2, "native_born", x=1, xpre=0, scen=9)

xhyp1.2 <- cfName(xhyp1.2, "female", scen=10) 
xhyp1.2 <-cfChange(xhyp1.2, "female", x=1, xpre=0, scen=10)

xhyp1.2 <- cfName(xhyp1.2, "age", scen=11) 
xhyp1.2 <-cfChange(xhyp1.2, "age", x=40, xpre=0, scen=11)

xhyp1.2 <- cfName(xhyp1.2, "education2", scen=12) 
xhyp1.2 <-cfChange(xhyp1.2, "education2", x=5, xpre=0, scen=12)

xhyp1.2 <- cfName(xhyp1.2, "highinc", scen=13) 
xhyp1.2 <-cfChange(xhyp1.2, "highinc", x=1, xpre=0, scen=13)

xhyp1.2 <- cfName(xhyp1.2, "medinc", scen=14) 
xhyp1.2 <-cfChange(xhyp1.2, "medinc", x=1, xpre=0, scen=14)

xhyp1.2 <- cfName(xhyp1.2, "black_caribbean", scen=15) 
xhyp1.2 <-cfChange(xhyp1.2, "black_caribbean", x=1, xpre=0, scen=15)

xhyp1.2 <- cfName(xhyp1.2, "indian", scen=16) 
xhyp1.2 <-cfChange(xhyp1.2, "indian", x=1, xpre=0, scen=16)

xhyp1.2 <- cfName(xhyp1.2, "pakistani", scen=17) 
xhyp1.2 <-cfChange(xhyp1.2, "pakistani", x=1, xpre=0, scen=17)

xhyp1.2 <- cfName(xhyp1.2, "bangladeshi", scen=18) 
xhyp1.2 <-cfChange(xhyp1.2, "bangladeshi", x=1, xpre=0, scen=18)

xhyp1.2 <- cfName(xhyp1.2, "vote_duty", scen=19) 
xhyp1.2 <-cfChange(xhyp1.2, "vote_duty", x=5, xpre=1, scen=19)

xhyp1.2 <- cfName(xhyp1.2, "efficacy", scen=20) 
xhyp1.2 <-cfChange(xhyp1.2, "efficacy", x=10, xpre=0, scen=20)

xhyp1.2 <- cfName(xhyp1.2, "democ_satis", scen=21) 
xhyp1.2 <-cfChange(xhyp1.2, "democ_satis", x=3, xpre=0, scen=21)

xhyp1.2 <- cfName(xhyp1.2, "trust_parliament", scen=22) 
xhyp1.2 <-cfChange(xhyp1.2, "trust_parliament", x=10, xpre=0, scen=22)

xhyp1.2
yhyp1.2<-logitsimfd(xhyp1.2, simbetas, ci=0.90)
yhyp1.2


##Graph
trace1 <- ropeladder(
  x=yhyp2.2$pe,
  lower=yhyp2.2$lower,
  upper=yhyp2.2$upper,
  labels=c("Societal Discrimination", "Political Discrimination",
  "Worship Attendance", "Political Interest", "Political Knowledge", "Party ID (Yes=1)", 
  "Identity (Brit=2)", "English (Main Lang)", "Native Born", "Female", 
  "Age", "Education", "High Income", "Med Income",
  "Black Caribbean", "Indian", "Pakistani", "Bangladeshi",
  "Vote Duty", "Political Efficacy","Democratic Satisfaction","Trust Parliament"),
  col= c("Black","Black","Black","Black","Black","Black","grey60","grey60","Black","Black","Black", "grey60", "grey60", "grey60", "grey60", "Black", "Black", "Black", "Black", "Black", "grey60", "Black"),
  pch=c(16,16,16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16),
  cex=2, 
  plot=1)

trace2 <- ropeladder(
  x=yhyp1.2$pe,
  lower=yhyp1.2$lower,
  upper=yhyp1.2$upper,
  col= c("Black","Black","Black","Black","Black","Black","grey60","Black","Black","Black","Black", "grey60", "grey60", "Black", "grey60", "Black", "Black", "Black", "Black", "Black", "grey60", "grey60"),
  pch=c(16,16,16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16),
  cex=2, 
  plot=2)

vertmark <- linesTile(x = c(0,0),
                      y = c(0,1),
                      lty = "dashed",
                      col= "Black",
                      plot = 1:2)

# Widen space between titles on left
trace1$entryheight <- .10
trace2$entryheight <- .10

at.x1 <- seq(-.4, .5,.1)
at.x2 <- seq(-.4, .5,.1)

xmin1<- -.4 
xmax1<-.5

xmin2<- -.4
xmax2<- .5

tc<-tile(trace1,
         trace2,
         vertmark,
         xaxis=list(at1 = at.x1, 
                    at2= at.x2, 
                    cex=.8),
         topaxis=list(at = at.x1, 
                      at2= at.x2, 
                      cex=.8,
                      add = rep(TRUE,2)),
         undertitle = list(labels="Change in Predicted Probability (Min-Max)", cex=.8, x=1),
         topaxistitle = list(labels=c("Local Election", "General Election"), cex=1),
         output = list(file = "EMBES_figure4", type="pdf", width=10), 
         frame=FALSE,
         gridlines=list(type="t", col="grey65")
)



###################################################################
##Table 3 (Impact of Street and Gov Disc on Electoral Particiation)
###################################################################

#Local
formula<-(voted2010_local ~ 
            street_racedisc + poldisc_gov_only+
            relatt_oth_r + pol_interest + polknowledge + partyid +
            identity+ english+ native_born+ female+ age+ education2+ highinc+ 
            medinc+ misinc+ black_caribbean+ indian+ pakistani+ bangladeshi+
            vote_duty+ efficacy+ democ_satis+ trust_parliament)


mdata<-extractdata(formula, data, na.rm=TRUE)

k <- length(names(mdata))
for (k in 1:k) {
  print(class(mdata[,k]))
}
for (k in 1:k) {
  mdata[,k] <- as.numeric(mdata[,k])
}

logit.result_local_robust<- glm(formula, family=binomial(link=logit), data=mdata)
summary(logit.result_local_robust)


#General
formula<-(voted2010 ~ 
            street_racedisc + poldisc_gov_only+
            relatt_oth_r + pol_interest + polknowledge + partyid +
            identity+ english+ native_born+ female+ age+ education2+ highinc+ 
            medinc+ misinc+ black_caribbean+ indian+ pakistani+ bangladeshi+
            vote_duty+ efficacy+ democ_satis+ trust_parliament)


mdata<-extractdata(formula, data, na.rm=TRUE)

k <- length(names(mdata))
for (k in 1:k) {
  print(class(mdata[,k]))
}
for (k in 1:k) {
  mdata[,k] <- as.numeric(mdata[,k])
}

logit.result_general_robust<- glm(formula, family=binomial(link=logit), data=mdata)
summary(logit.result_general_robust)


modellabels<-c( "Street Disc (Societal)", "Government Disc (Political)",
               "Worship Attendance", "Political Interest", "Political Knowledge", "Party ID (Yes=1)", 
               "Identity (Brit=2)", "English (Main Lang)", "Native Born", "Female", 
               "Age", "Education", "High Income", "Med Income", "Missing Income",
               "Black Caribbean", "Indian", "Pakistani", "Bangladeshi",
               "Vote Duty", "Political Efficacy","Democratic Satisfaction","Trust Parliament")


dvlabel1<-c("Local Election", "General Election")


stargazer(logit.result_local_robust, logit.result_general_robust, style="APSR", covariate.labels = modellabels, out.header=T,
          model.numbers = TRUE, 
          title = "",
          dep.var.labels = dvlabel1)




###############################################################
##Table 5: Impact of Discrimination on Ethnic-Based Engagement
###############################################################

#Full Model
formula<-(ethnic_active~ socdisc + poldisc +
            relatt_oth_r + pol_interest+ polknowledge+ partyid + identity+ english+ native_born+
            female+ age+ education2+ highinc+ medinc+ misinc+ black_caribbean+ indian+ pakistani+ bangladeshi + efficacy+ democ_satis+ trust_parliament)


mdata<-extractdata(formula, data, na.rm=TRUE)

k <- length(names(mdata))
for (k in 1:k) {
  print(class(mdata[,k]))
}
for (k in 1:k) {
  mdata[,k] <- as.numeric(mdata[,k])
}

logit.result3<- glm(formula, family=binomial(link=logit), data=mdata)
summary(logit.result3)

pe.3 <- logit.result3$coefficients  
vc.3 <- vcov(logit.result3)


# simulation
sims <- 10000
simbetas <- mvrnorm(sims, pe.3, vc.3)

xhyp3 <- cfMake(formula, mdata, nscen=2, f=mean)

xhyp3 <- cfName(xhyp3, "socdisc", scen=1) 
xhyp3 <-cfChange(xhyp3, "socdisc", x=9, xpre=0, scen=1)

xhyp3 <- cfName(xhyp3, "poldisc", scen=2) 
xhyp3 <-cfChange(xhyp3, "poldisc", x=12, xpre=0, scen=2)

xhyp3
yhyp3<-logitsimfd(xhyp3, simbetas, ci=0.90)
yhyp3


#Model Without Efficacy, Gov. Trust, and Dem. Satisfaction
formula<-(ethnic_active~ socdisc + poldisc +
            relatt_oth_r + pol_interest+ polknowledge+ partyid + identity+ english+ native_born+
            female+ age+ education2+ highinc+ medinc+ misinc+ black_caribbean+ indian+ pakistani+ bangladeshi)


mdata<-extractdata(formula, data, na.rm=TRUE)

k <- length(names(mdata))
for (k in 1:k) {
  print(class(mdata[,k]))
}
for (k in 1:k) {
  mdata[,k] <- as.numeric(mdata[,k])
}

logit.result4<- glm(formula, family=binomial(link=logit), data=mdata)
summary(logit.result4)

modellabels2<-c ("Societal Disc", "Political Disc",
                 "Worship Attendance", "Political Interest", "Political Knowledge", "Party ID (Yes=1)", 
                 "Identity (Brit=2)", "English (Main Lang)", "Native Born", "Female", 
                 "Age", "Education", "High Income", "Med Income", "Missing Income",
                 "Black Caribbean", "Indian", "Pakistani", "Bangladeshi", "Political Efficacy", "Democratic Satisfaction", "Trust Parliament")


dvlabel3<-c("Ethnic-Based Participation")


stargazer(logit.result3,logit.result4, style="APSR", covariate.labels = modellabels2, out.header=T,
          model.numbers = TRUE, 
          title = "",
          dep.var.labels = dvlabel3)



##################################################################
##Table 6: Impact of Street & Gov Disc on Ethnic-Based Engagement
##################################################################
formula<-(ethnic_active~ street_racedisc + poldisc_gov_only+
            relatt_oth_r + pol_interest+ polknowledge+ partyid + identity+ english+ native_born+
            female+ age+ education2+ highinc+ medinc+ misinc+ black_caribbean+ indian+ pakistani+ bangladeshi + efficacy+ democ_satis+ trust_parliament)

mdata<-extractdata(formula, data, na.rm=TRUE)

k <- length(names(mdata))
for (k in 1:k) {
  print(class(mdata[,k]))
}
for (k in 1:k) {
  mdata[,k] <- as.numeric(mdata[,k])
}

logit.result5<- glm(formula, family=binomial(link=logit), data=mdata)
summary(logit.result5)

modellabels3<-c ("Street Disc (Societal)", "Government Disc (Political)",
                 "Worship Attendance", "Political Interest", "Political Knowledge", "Party ID (Yes=1)", 
                 "Identity (Brit=2)", "English (Main Lang)", "Native Born", "Female", 
                 "Age", "Education", "High Income", "Med Income", "Missing Income",
                 "Black Caribbean", "Indian", "Pakistani", "Bangladeshi", "Political Efficacy", "Democratic Satisfaction", "Trust Parliament")


dvlabel4<-c("Ethnic-Based Participation")


stargazer(logit.result5, style="APSR", covariate.labels = modellabels3, out.header=T,
          model.numbers = TRUE, 
          title = "",
          dep.var.labels = dvlabel4)


##################################################
##Table 4: The Impact of Disc on Identity Choice
#################################################
library(nnet)
library(broom)
library(effects)

formula<- (identity ~ socdisc + poldisc +
           relatt_oth_r + pol_interest+ polknowledge+ partyid+ english+ native_born+
             female+ age+ education2+ highinc+ medinc+ misinc+ black_caribbean+ indian+ pakistani+ bangladeshi + efficacy+ democ_satis+ trust_parliament)

mdata <- extractdata(formula, data, na.rm=TRUE) 

k <- length(names(mdata))
for (k in 1:k) {
  print(class(mdata[,k]))
}
for (k in 1:k) {
  mdata[,k] <- as.numeric(mdata[,k])
}

mdata$identity <- factor(mdata$identity, 
                         levels= c(0, 1, 2),
                             labels=c("BlackAsian", "Both", "British"))


mlogitresult <- multinom(identity ~ socdisc + poldisc +
                         relatt_oth_r + pol_interest+ polknowledge+ partyid + english+ native_born+
                           female+ age+ education2+ highinc+ medinc+ misinc+ black_caribbean+ indian+ pakistani+ 
                           bangladeshi + efficacy+ democ_satis+ trust_parliament, Hess = TRUE, data=mdata)

summary(mlogitresult)


#Models without Efficacy, Dem. Satis, and GOV. Trust
formula<- (identity ~ socdisc + poldisc +
             relatt_oth_r + pol_interest+ polknowledge+ partyid + english+ native_born+
             female+ age+ education2+ highinc+ medinc+ misinc+ black_caribbean+ indian+ pakistani+ bangladeshi)

mdata <- extractdata(formula, data, na.rm=TRUE) 

k <- length(names(mdata))
for (k in 1:k) {
  print(class(mdata[,k]))
}
for (k in 1:k) {
  mdata[,k] <- as.numeric(mdata[,k])
}

mdata$identity <- factor(mdata$identity, 
                         levels= c(0, 1, 2),
                         labels=c("BlackAsian", "Both", "British"))


mlogitresult2 <- multinom(identity ~ socdisc + poldisc +
                           relatt_oth_r + pol_interest+ polknowledge+ partyid + english+ native_born+
                           female+ age+ education2+ highinc+ medinc+ misinc+ black_caribbean+ indian+ pakistani+ 
                           bangladeshi, Hess = TRUE, data=mdata)

summary(mlogitresult2)


#Table w/ Stargazer
tab_names <-c("Societal Disc","Political Disc", "Worship Attendance", "Political Interest", "Political Knowledge", "Party ID (Yes=1)", 
               "English (Main Lang)", "Native Born", "Female", 
              "Age", "Education", "High Income", "Med Income", "Missing Income",
              "Black Caribbean", "Indian", "Pakistani", "Bangladeshi", "Political Efficacy", "Democratic Satisfaction", "Trust Parliament")

covariate.labels = tab_names

sg_out <- stargazer(mlogitresult, mlogitresult2, covariate.labels = tab_names, style="apsr",title = "Identity")

write(sg_out, file=paste("EMBES_Table4.tex",sep="/"))


########################################################################################
##Figure 5: Relationship btwn Disc and Identity (Imported First Differences From Stata Using Margins Function)
########################################################################################

#Identifying Asian/Black relative to British 
yhype5p <- c(0.3065, 0.2568) 
yhype5l<- c(0.1642, 0.0653)
yhype5u<- c(0.4487, 0.4483)

#Identifying Both Relative to Asian/Black
yhype5.2p <- c(-0.2611, -0.1465) 
yhype5.2l<- c(-0.3829, -0.3144)
yhype5.2u<- c(-0.1392, -0.005)

##Graphing 
trace1 <- ropeladder(
  x=yhype5p,
  lower=yhype5l,
  upper=yhype5u,
  labels=c("Societal Discrimination", "Political Discrimination"), 
  col="Black", 
  pch=c(16,16),
  cex=1.8, 
  plot=1)

trace2 <- ropeladder(
  x=yhype5.2p,
  lower=yhype5.2l,
  upper=yhype5.2u,
  col="Black",
  pch=c(16,16),
  cex=1.8,
  plot=2)

vertmark <- linesTile(x = c(0,0),
                      y = c(0,1),
                      lty = "dashed",
                      col= "Black",
                      plot = 1:2)

trace1$entryheight <- .25
trace2$entryheight <- .25

at.x1 <- seq(-.1,.4,.1)
at.x2 <-seq(-.4,.1,.1)

tc<-tile(trace1,
         trace2,
         vertmark,
         xaxis=list(at1 = at.x1, 
            at2= at.x2, 
         cex=.8), 
        topaxis=list(at1 = at.x1, 
                     at2= at.x2, 
                     cex=.8,
                     add = rep(TRUE,2)),
         undertitle = list(labels="Change in Predicted Probability (Min-Max)", cex=.8, x=1),
         plottitle= list(labels1= "Black/Asian vs. British", 
         labels2= "Equally Both vs. Black/Asian"),
         output = list(file = "EMBES_figure5", width=12, type="pdf"), 
         frame=FALSE,
         gridlines=list(type="x", col="grey65")
)




##############################################################
##Figure 6: Relationship btwn Disc and Ethnic-Based Engagmenet
##############################################################
formula<-(ethnic_active~ socdisc + poldisc +
            relatt_oth_r + pol_interest+ polknowledge+ partyid + identity+ english+ native_born+
            female+ age+ education2+ highinc+ medinc+ misinc+ black_caribbean+ indian+ pakistani+ bangladeshi + efficacy+ democ_satis+ trust_parliament)


mdata<-extractdata(formula, data, na.rm=TRUE)

k <- length(names(mdata))
for (k in 1:k) {
  print(class(mdata[,k]))
}
for (k in 1:k) {
  mdata[,k] <- as.numeric(mdata[,k])
}

logit.result3<- glm(formula, family=binomial(link=logit), data=mdata)
summary(logit.result3)

pe.3.2 <- logit.result3$coefficients  
vc.3.2 <- vcov(logit.result3)

### simulation
sims <- 10000
simbetas <- mvrnorm(sims, pe.3.2, vc.3.2)

xhyp3.2 <- cfMake(formula, mdata, nscen=21, f=mean)

xhyp3.2 <- cfName(xhyp3.2, "socdisc", scen=1) 
xhyp3.2 <-cfChange(xhyp3.2, "socdisc", x=9, xpre=0, scen=1)

xhyp3.2 <- cfName(xhyp3.2, "poldisc", scen=2) 
xhyp3.2 <-cfChange(xhyp3.2, "poldisc", x=12, xpre=0, scen=2)

xhyp3.2 <- cfName(xhyp3.2, "relatt_oth_r", scen=3) 
xhyp3.2 <-cfChange(xhyp3.2, "relatt_oth_r", x=5, xpre=0, scen=3)

xhyp3.2 <- cfName(xhyp3.2, "pol_interest", scen=4) 
xhyp3.2 <-cfChange(xhyp3.2, "pol_interest", x=4, xpre=0, scen=4)

xhyp3.2 <- cfName(xhyp3.2, "polknowledge", scen=5) 
xhyp3.2 <-cfChange(xhyp3.2, "polknowledge", x=5, xpre=0, scen=5)

xhyp3.2 <- cfName(xhyp3.2, "partyid", scen=6) 
xhyp3.2 <-cfChange(xhyp3.2, "partyid", x=1, xpre=0, scen=6)

xhyp3.2 <- cfName(xhyp3.2, "identity", scen=7) 
xhyp3.2 <-cfChange(xhyp3.2, "identity", x=2, xpre=0, scen=7)

xhyp3.2 <- cfName(xhyp3.2, "english", scen=8) 
xhyp3.2 <-cfChange(xhyp3.2, "english", x=1, xpre=0, scen=8)

xhyp3.2 <- cfName(xhyp3.2, "native_born", scen=9) 
xhyp3.2 <-cfChange(xhyp3.2, "native_born", x=1, xpre=0, scen=9)

xhyp3.2 <- cfName(xhyp3.2, "female", scen=10) 
xhyp3.2 <-cfChange(xhyp3.2, "female", x=1, xpre=0, scen=10)

xhyp3.2 <- cfName(xhyp3.2, "age", scen=11) 
xhyp3.2 <-cfChange(xhyp3.2, "age", x=40, xpre=18, scen=11)

xhyp3.2 <- cfName(xhyp3.2, "education2", scen=12) 
xhyp3.2 <-cfChange(xhyp3.2, "education2", x=5, xpre=0, scen=12)

xhyp3.2 <- cfName(xhyp3.2, "highinc", scen=13) 
xhyp3.2 <-cfChange(xhyp3.2, "highinc", x=1, xpre=0, scen=13)

xhyp3.2 <- cfName(xhyp3.2, "medinc", scen=14) 
xhyp3.2 <-cfChange(xhyp3.2, "medinc", x=1, xpre=0, scen=14)

xhyp3.2 <- cfName(xhyp3.2, "black_caribbean", scen=15) 
xhyp3.2 <-cfChange(xhyp3.2, "black_caribbean", x=1, xpre=0, scen=15)

xhyp3.2 <- cfName(xhyp3.2, "indian", scen=16) 
xhyp3.2 <-cfChange(xhyp3.2, "indian", x=1, xpre=0, scen=16)

xhyp3.2 <- cfName(xhyp3.2, "pakistani", scen=17) 
xhyp3.2 <-cfChange(xhyp3.2, "pakistani", x=1, xpre=0, scen=17)

xhyp3.2 <- cfName(xhyp3.2, "bangladeshi", scen=18) 
xhyp3.2 <-cfChange(xhyp3.2, "bangladeshi", x=1, xpre=0, scen=18)

xhyp3.2 <- cfName(xhyp3.2, "efficacy", scen=19) 
xhyp3.2 <-cfChange(xhyp3.2, "efficacy", x=10, xpre=0, scen=19)

xhyp3.2 <- cfName(xhyp3.2, "democ_satis", scen=20) 
xhyp3.2 <-cfChange(xhyp3.2, "democ_satis", x=3, xpre=0, scen=20)

xhyp3.2 <- cfName(xhyp3.2, "trust_parliament", scen=21) 
xhyp3.2 <-cfChange(xhyp3.2, "trust_parliament", x=10, xpre=0, scen=21)

xhyp3.2
yhyp3.2<-logitsimfd(xhyp3.2, simbetas, ci=0.90)
yhyp3.2


trace1 <- ropeladder(
  x=yhyp3.2$pe,
  lower=yhyp3.2$lower,
  upper=yhyp3.2$upper,
  labels=c("Societal Discrimination", "Political Discrimination",
           "Worship Attendance", "Political Interest", "Political Knowledge", "Party ID (Yes=1)", 
           "Identity (Brit=2)", "English (Main Lang)", "Native Born", "Female", 
           "Age", "Education", "High Income", "Med Income",
           "Black Caribbean", "Indian", "Pakistani", "Bangladeshi", 
           "Political Efficacy","Democratic Satisfaction","Trust Parliament"),
  col= c("Black","Black","Black","Black","Black","Black","Black","grey60","Black","grey60","grey60", "Black", "grey60", "grey60", "grey60", "Black", "Black", "grey60", "Black", "grey60", "grey60"),
  pch=c(16,16,16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16),
  cex=2)


vertmark <- linesTile(x = c(0,0),
                      y = c(0,1),
                      lty = "dashed",
                      col= "Black",
                      plot = 1)

trace1$entryheight <- .10

at.x1 <- seq(-.2, .4,.1)

xmin1<- -.2 
xmax1<-.4

xmin2<- -.2
xmax2<- .4

tc<-tile(trace1,
         vertmark,
         xaxis=list(at1 = at.x1, 
                    cex=.8),
         topaxis=list(at = at.x1, 
                      cex=.8,
                      add = rep(TRUE,1)),
         undertitle = list(labels="Change in Predicted Probability (Min-Max)", cex=.8, x=.5),
        topaxistitle = list(labels=c("Ethnic-Based Engagement"), cex=1),
         output = list(file = "EMBES_figure6", type="pdf", width=8), 
         frame=FALSE,
         gridlines=list(type="t", col="grey65")
)

